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Motivated by the properties of the iron chalcogenides, we study the phase diagram of a generalized Heisen¬ 
berg model with frustrated bilinear-biquadratic interactions on a square lattice. We identify zero-temperature 
phases with antiferroquadrupolar and Ising-nematic orders. The effects of quantum fluctuations and interlayer 
couplings are analyzed. We propose the Ising-nematic order as underlying the structural phase transition ob¬ 
served in the normal state of FeSe, and discuss the role of the Goldstone modes of the antiferroquadrupolar order 
for the dipolar magnetic fluctuations in this system. Our results provide a considerably broadened perspective 
on the overall magnetic phase diagram of the iron chalcogenides and pnictides, and are amenable to tests by 
new experiments. 
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Introduction. Because superconductivity develops near 
magnetic order in most of the iron pnictides and chalco¬ 
genides, it is important to understand the nature of their 
magnetism. The iron pnictide families typically have parent 
compounds that show a collinear ( 7 r, 0) antiferromagnetic or¬ 
der m. Lowering the temperature in the parent compounds 
gives rise to a tetragonal-to-orthorhombic distortion, and the 
temperature T s for this structural transition is either equal to 
or larger than the Neel transition temperature 7’y. A likely ex¬ 
planation for T s is an Ising-nematic transition at the electronic 
level. It was recognized from the beginning that models with 
quasi-local moments and their frustrated Heisenberg ,1\ — ./ 2 
interactions 13 feature such an Ising-nematic transition SI). 
Similar conclusions have subsequently been reached in mod¬ 
els that are based on Fermi-surface instabilities (7). 

The magnetic origin for the nematicity fits well with the 
experimental observations of the spin excitation spectrum ob¬ 
served in the iron pnictides. Inelastic neutron scattering ex¬ 
periments ED in the parent iron pnictides have revealed a 
low-energy spin spectrum whose equal-intensity contours in 
the wavevector space form ellipses near (±7r, 0) and (0, ±7r). 
At high energies, spin-wave-like excitations are observed all 
the way to the boundaries of the underlying antiferromag¬ 
netic Brillouin zone 0. These features are well captured by 
Heisenberg models with the frustrated Ji — J 2 interactions 
and biquadratic couplings mm, although at a refined level 
it is also important to incorporate the damping provided by the 
coherent itinerant fermions near the Fermi energy ED. 

Experiments in bulk FeSe do not seem to fit into this frame¬ 
work. FeSe is one of the canonical iron chalcogenide su¬ 
perconductors Hi. In the single-layer limit, it currently 
holds particular promise towards a very high T c superconduc¬ 
tivity QME1 driven by strong correlations ini- In the bulk 
form, this compound displays a tetragonal-to-orthorhombic 
structural transition, with T s « 90 K, but no Neel transi¬ 
tion has been detected fl8|j2TI . This distinction has been in¬ 


terpreted as pointing towards the failure of the magnetism- 
based origin for the structural phase transition |20][2Tj|. At the 
same time, experiments have also revealed another aspect of 
the emerging puzzle. The structural transition clearly induces 
dipolar magnetic fluctuations |[20l [211. 

In this Letter, we show that a generalized Heisenberg model 
with frustrated bilinear-biquadratic couplings on a square lat¬ 
tice contains a phase with both a [it, 0) antiferroquadrupolar 
(AFQ) order and an Ising-nematic order. The model in this 
regime displays a finite-temperature transition into the Ising- 
nematic order and, in the presence of interlayer couplings, 
also a finite-temperature transition into the AFQ order. We 
suggest that such physics is compatible with the aforemen¬ 
tioned and related properties of FeSe. The Goldstone modes 
of the AFQ order are responsible for the onset of dipolar mag¬ 
netic fluctuations near the wave-vector ( 7 r, 0), which is exper¬ 
imentally testable. 

Model. We consider a spin Hamiltonian with S 1 on a 
two-dimensional (2D) square lattice: 

h = \J2 { • s 7 + K ^i ■ Sj) 2 } , a) 

i,6 n 

where j = i + S n , and S n connects site i and its n-th near¬ 
est neighbor sites with n = 1, 2, 3. Here ./„ and K n are, 
respectively, the bilinear and biquadratic couplings between 
the n-th nearest neighbor spins. For iron pnictides and iron 
chalcogenides, the local moments describe the spin degrees 
of freedom associated with the incoherent part of the elec¬ 
tronic excitations and reflect the bad-metal behavior of these 
systems U] HI $-lbJ. A nonzero J 3 is believed to be impor¬ 
tant for the iron chalcogenides, especially FeTe E3 . The bi¬ 
quadratic couplings K n are expected to play a significant role 
in multi-orbital systems with moderate Hund’s coupling (23). 
The nearest neighbor coupling K 1 was included in previous 
studies mm to understand the strong anisotropic spin ex¬ 
citations in the Ising-nematic ordered phase, where the ground 
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FIG. 1. Momentum distribution of the dipolar (top row) and 
quadrupolar (bottom row) magnetic structure factors at A 2 = — 1 
(in (a) and (c)) and A '2 = 1.5 (in (b) and (d)), respectively. Here, 
Ji = J 2 = 1, J 3 = K 3 = 0, and A'i = — 1. The calculations 
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are done on a 40x40 lattice at T/|Ai|=0.01 with up to 10 s Monte 
Carlo steps. In (d), besides the leading AFQ correlations at ( 7 r, 0) 
and (0, 7 r), subleading FQ correlations are observed at finite temper¬ 
atures; as the temperature is lowered, the former is enhanced whereas 
the latter is is diminished. 


state has a ( 7 r, 0)/(0, ir) antiferromagnetic (AFM) order. With 
the goal of searching for the new physics that could describe 
the properties of FeSe, in this work, we take these couplings as 
variables and study the pertinent unusual phases in the phase 
diagram. In the following, to simplify the discussion on the 
relevant AFM and AFQ phases, we take K t = — 1 and use 
|A'i| as the energy unit. Note that a moderately positive A'i 
in the presence of further-neighbor K n couplings will lead to 
similar results, but a A'i coupling alone in the absence of the 
latter will not generate the physics discussed below. 

Some general considerations are in order. For S ^ 1 

(Si • S,.) 2 = • Q j - iSi • S, + *S 2 S 2 , (2) 


where Q, is a quadrupolar operator with five components 

Qf- 3 * 2 = ^ra 2 +(^) 2 -2 (stnQf- y2 = (Af) 2 - 
( S y) 2 , or* = Sfsf + s?sf, QT = svsf + s*si 


and Q\ x = S\S X + SfSf. Therefore, the biquadratic 
interaction may promote a long-range ferroquadrupolar- 
antiferroquadrupolar (FQ-AFQ) order. With the aforemen¬ 
tioned motivation, we are interested in a ( 7 r, 0) AFQ order, 
which would break the C 4 symmetry and should yield an 
Ising-nematic order parameter. 

Low-temperature phase diagram of the classical spin 
model. We first study the model in Eq. |T]) for classical spins. 
For simplicity, we discuss the case A 3 = 0. We have calcu¬ 
lated the dipolar and quadrupolar magnetic structure factors 
via Monte Carlo simulations using the standard Metropolis 


FIG. 2. Low-temperature phase diagrams of the classical bilinear- 
biquadratic Heisenberg model at (a): Ji = J 2 , J 3 = A 3 = 0 and 
(b): Ji = A 3 = 0, J 2 = J 3 . Both are shown at T/|Ai| = 0.01. 
Dashed lines show finite-temperature crossovers between different 
orders. The dominant order in each regime is labeled. In each case, 
the solid line shows the mean-field phase boundary at T = 0. 


algorithm. (24]] Representative results for the structure factor 
data are shown in Fig. [T] for J 3 = 0 and J\ = J 2 . The 
two cases, corresponding to different values of A' 2 , show, re¬ 
spectively, dominant ferroquadrupolar (FQ) and (n, 0) AFQ 
correlations, for the finite-size systems studied and at a very 
low temperature T/|A'i| =0.01. 

Overall, as shown in Fig. |2|a), we find that there are large 
regimes in the phase diagram in which the FQ and ( 7 r, 0) AFQ 
moments are almost ordered, while the dipolar moments co¬ 
existing with the FQ/AFQ moments are very weakly corre¬ 
lated. Hence in these regimes, the dominant low-temperature 
order is the FQ/AFQ one. In between these, there is a regime 
in which the dominant correlation occurs in the (n, 0) AFM 
channel. 

Similar results for the case of Ji =0 and J 2 = J 3 are 
shown in Fig. [2jb). A large regime with dominating FQ or 
( 7 r, 0) AFQ correlations is also found. The difference from 
the case of J 3 = 0 and J\ = J 2 occurs in the regime with 
dominant AFM correlations, for which the wavevector is now 
(±7t/2, ±7t/2) as relevant to the FeTe compound. 

For 2D systems, thermal fluctuations will ultimately (in the 
thermodynamic limit) destroy any order that breaks a con¬ 
tinuous global symmetry at any nonzero temperature ll25l . 
The dashed lines in Fig.[2]therefore mark crossovers between 
regimes with different dominant correlations. At T = 0, on 
the other hand, genuine FQ/AFQ can occur in our model on 
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FIG. 3. Temperature dependence of the Ising-nematic order param¬ 
eters <ti and (72 at (a) JT = J 2 = J 3 = 0, AT = — I, and AT = 1 
and (b) Ji = 0, J 2 = J 3 = 0.5, AT = —1. and AT = 2. In both 
cases the dominant part of the Ising-nematic order is a 2 , which is 
associated with the AFQ order. 


the square lattice. We have therefore also analyzed the mean- 
field phase diagrams at T = 0. The resulting phase boundary 
is shown in each case as a solid line in Fig. [2] The results are 
compatible with the crossovers identified at low but nonzero 
temperatures. For the case of J 3 = 0 and J\ = ■!■>, shown in 
Fig.[2|a), the phase on the left of the solid line has a mixture of 
an AFM phase ordered at q = ( 7 r, 0)/(0, 7 r) and a FQ phase. 
The phase on the right of the solid line has an AFQ phase or¬ 
dered at q = ( 77 , 0)/(0,7r). Note that in the classical limit, 
the spins are treated as 0(3) vectors, and should always be 
ordered at zero temperature. We find that in the AFQ phase, 
the spins can be ordered at a wavevector ( q , 7r)/(7r, q) for ar¬ 
bitrary q, with an infinite degeneracy. (26) Such a frustration 
would likely stabilize a purely AFQ ground state when quan¬ 
tum fluctuations are taken into account (see below). For the 
case of J\ =0 and J 2 = J 3 , shown in Fig. |2fb), the mean- 
field result also yields FQ or (n, 0) AFQ, respectively, to the 
left or right of the solid line. However, the wave vector for 
the AFM orders that mix, respectively, with the FQ and (tt, 0) 
AFQ order has become (± 7 t/ 2 , ±7r/2). |[26i 

Similar to the ( 7 r, 0 ) AFM state, the ( 7 r, 0) AFQ phase 
breaks the lattice C 4 rotational symmetry. An accompanying 
Ising-nematic transition is to be expected, and should develop 
at nonzero temperatures even in two dimensions. We define 
the general Ising-nematic operators as follows: 

= (Si • S i+a ) n - (Si • S i+ #) n , (3) 

where n = 1,2. We also introduce the quadrupolar Qa/b to 
be the linear superposition of Q(7t, 0)/(0, it), with the ratios 
of their coefficients to be ±1 respectively. From Eq. ([2]), we 
see that for quantum spins, the Ising-nematic order associated 
with Q should be seen in both 04 and cr 2 . For classical spins, 
since Q, • Q j = 2 (S; ■ Sj) 2 — |S 2 S 2 , only er 2 will manifest 
Qt • Q /i • This allows us to determine the origin of the Ising- 
nematic order in the AFQ+AFM phase. As shown in Fig.|3ja), 



FIG. 4. Calculated spin excitation spectrum in the ( 7 r, 0) AFQ phase 
of the quantum 5 = 1 model. We have taken J\ = J 2 = 0.25, 
J 3 = 0, K 2 = 0.5, and A 3 = —0.3. The color codes the dynamical 
spin dipolar structure factor, yj 5| “(q, lj). 


for the Ai — K 2 model, a 2 is ordered at T a /\K\\ ~ 0.38 but 
o\ ~ 0 for any T > 0. Likewise from Fig. [3j b), in the case 
Ji = 0 and J ‘2 = J 3 , the dominant Ising nematic order pa¬ 
rameter is 02 for T < T a sa 0.9, and oq never becomes sub¬ 
stantial down to the lowest temperature T = 10 -4 accessible 
to our numerical simulation. These indicate that the Ising- 
nematic order in the AFQ+AFM phase is associated with the 
anisotropic spin quadrupolar fluctuations. 

The quantum spin models. The AFQ phase and the as¬ 
sociated Ising-nematic transition are features of the general¬ 
ized J — K model for both classical and quantum spins. To 
consider the effect of quantum fluctuations, we consider the 
case of S = 1. We study its ground-state properties via a 
semiclassical variational approach by using an SU(3) repre¬ 
sentation G2, and identify parameter regimes that stabilize 
the AFQ phase. We further study the spin excitations in the 
AFQ phase with the ordering wavevector qq = (714 0) using 
a flavor-wave theory. lf26ll Because the AFQ order breaks the 
continuous spin-rotational invariance, the Goldstone modes 
will have a nonzero dipolar matrix element Ell [28). To 
explicitly demonstrate this, we calculate the dynamical spin 
dipolar structure factor Sfffq. ui) near qq, which is shown 
in Fig. [4] Therefore, the development of the AFQ order is ac¬ 
companied by a sharp rise in the dynamical spin dipolar corre¬ 
lations centered around the wavevector ( 7 r, 0) (and symmetry- 
related wave vectors). 

Coupling to itinerant fermions and interaction between lay¬ 
ers. One additional effect of the quantum fluctuations is that 
it can suppress the weak AFM order when the dominant or¬ 
der is AFQ. We discuss one source of such an effect, which 
is the coupling of the order parameters to the coherent itiner¬ 
ant fermions. The effect of coupling to the itinerant fermions 
can be treated as in Ref. (6| within an effective Ginzburg- 
Landau action, and is briefly discussed in the Supplemen¬ 
tal Material l26l . When only the (7r,0) AFM order and the 
Ising-nematic order are present, the coupling to the itinerant 
fermions will suppress the AFM and Ising-nematic order con¬ 
currently |29l . However, when the dominant order is AFQ, 
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FIG. 5. Skecthed phase diagram in terms of T and Ji/Ki- The 
dominant order may be either AFQ or AFM. The thinner dashed 
lines show the associated ordering temperautures Tafq and Tafm- 
A first-order transition (thicker dashed line) takes place at an in¬ 
termediate J 2 //T 2 coupling when the two transitions meet. The 
Ising-nematic transition (solid line) takes place at T a . There can 
be either a first-order Ising-nematic and AFM(AFQ) transition at 
T a = Tafm/afq , or two separate transitions. 


the coupling to the itinerant fermions can suppress the AFM 
order while retaining the stronger AFQ order and the associ¬ 
ated Ising-nematic order. 

When interlayer bilinear-biquadratic couplings are taken 
into account, a phase with a pure AFQ order can be stabi¬ 
lized at finite temperature. We can then discuss the evolution 
of the Ising-nematic transition as a function of the J 2 IK 2 ra¬ 
tio. Consider the case when a dominating J 2 stabilizes a (7r, 0) 
AFM order, which is accompanied by the Ising-nematic order 
parameter a\. For sufficiently large K- 2 , the AFQ order be¬ 
comes the dominant order, and the Ising-nematic order is pre¬ 
dominantly given by <72 ■ The schematic evolution between the 
two limits is illustrated in Fig. [5] We have illustrated the case 
with the quantum fluctuations having suppressed the weaker 
order. 

We stress that, such an evolution of the Ising-nematic tran¬ 
sition already occurs in the purely 2D model. Results from 
explicit calculations on the evolution of the transition temper¬ 
ature T a are shown in the Supplemental Material. |26| In the 
case of the Ising-nematic transition associated with a (ir, 0) 
AFM order, the interlayer couplings give rise to a nonzero 
Tafm < T„ (Refs.S®. Similarly, when the dominant or¬ 
der is a ( 7 r, 0) AFQ order, such couplings lead to a nonzero 
Tafq < T a . 

Implications for FeSe. General considerations suggest that 
the cases of spin 1 or spin 2 are pertinent to the iron-based 
materials EL Judging from the measured total spin spectral 
weight ID, the spin 1 case would be closer to the iron pnic- 
tides while the spin 2 case would be more appropriate for the 
iron chalcogenides. 

Accordingly, it is natural to propose that the normal state 
of FeSe realizes the phase whose ground state has the (7r, 0) 


AFQ order accompanied by the Ising-nematic order. In this 
picture, the structural transition at T s ~ 90 K corresponds 
to the concurrent Ising-nematic and AFQ transition, as illus¬ 
trated in Fig. [5] This picture explains why the structural phase 
transition is not accompanied by any static AFM order. At 
the same time, as soon as the AFQ order is developed, its 
Goldstone modes will contribute towards low-energy dipo¬ 
lar magnetic fluctuations. This is consistent with the onset of 
low-energy spin fluctuations observed in the NMR measure¬ 
ments Gol ED. It will clearly be important to explore such 
spin fluctuations using inelastic neutron scattering measure¬ 
ments. A quantitative comparison between the measured and 
calculated spin excitation spectra would allow estimates of the 
coupling constants J n and K n . The Goldstone modes may 
also be probed by magnetoresistance, and unusual features in 
this property have recently been reported lf30l . Finally, the 
Ising-nematic order is linearly coupled not only to the struc¬ 
tural anisotropy, but also to the orbital order. Similarly as for 
the iron pnictides ED, this would result in, for instance, the 
lifting of the d xz /d yz orbital degeneracy at the structural phase 
transition 13214341 . 

The phase diagrams given in Fig. [2] show that the AFQ 
region can be tuned to an AFM region. The nature of the 
AFM phase depends on the bilinear couplings. For a range 
of bilinear couplings, the nearby AFM phase has the order¬ 
ing wavevector (7r/2, 7 t/2). This provides a means to connect 
the magnetism of FeSe and FeTe [35]|36 j, which is of consid¬ 
erable interest to the on-going experimental efforts in study¬ 
ing the magnetism of the Se-doped FeTe series EL It also 
makes it natural to understand the development of magnetic 
order that seems to occur when FeSe is placed under a pres¬ 
sure on the order of 1 GPa 13814401 . Finally, we note that our 
results will serve as the basis to shed new light on the nematic 
correlations in the superconducting state ED®. 

Broader context. It is widely believed that understanding 
the magnetism in the iron chalcogenide FeTe, where the order¬ 
ing wavvector (n/2, it /2) has no connection with any Fermi- 
surface-nesting features ll35l [36), requires a local-moment 
picture. The proposal advanced here not only provides an 
understanding of the emerging puzzle on the magnetism in 
FeSe, but also achieves a level of commonality in the un¬ 
derlying microscopic interactions across these iron chalco¬ 
genides. Furthermore, the connection between the AFQ or¬ 
der and the (7r, 0) AFM order suggests that the local-moment 
physics, augmented by a coupling to the coherent itinerant 
fermions near the Fermi energy, places the magnetism of a 
wide range of iron-based superconductors in a unified frame¬ 
work. Since local-moment physics in bad metals reflects 
a proximity to correlation-induced electron localization, this 
unified perspective also signifies the importance of electron 
correlations 0 1441448) to the iron-based superconductors. 

Conclusions. To summarize, we have studied a general¬ 
ized Heisenberg model with frustrated bilinear-biquadratic 
interactions on a square lattice and find that the zero- 
temperature phase diagram stabilizes an antiferroquadrupo- 
lar order. The anisotropic spin quadrupolar fluctuations give 
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rise to a finite-temperature Ising-nematic transition. We pro¬ 
pose that the structural phase transition in FeSe corresponds 
to this Ising-nematic transition and is accompanied by an an- 
tiferroquadrupolar ordering. We suggest that inelastic neutron 
scattering experiments be carried out to explore the proposed 
Goldstone modes associated with the antiferroquadrupolar or¬ 
der. Our results provide a natural understanding for an emerg¬ 
ing puzzle on FeSe. More generally, the extended phase di¬ 
agrams advanced here considerably broaden the perspective 
on the magnetism and electron correlations of the iron-based 
superconductors. 

Note added. Recently, a study appeared that also em¬ 
phasized the local-moment-based magnetic physics for FeSe, 
but invoked a different mechanism based on a possible para¬ 
magnetic Ising-nematic ground state caused by J\-Ji frus¬ 
tration 1491 . A distinction of the mechanism advanced here 
is that the AFQ order yields Goldstone modes and therefore 
causes the onset of low-energy dipolar magnetic fluctuations. 
In addition, results from inelastic neutron scattering experi¬ 
ments in FeSe have appeared 15011511 . which verify the (7r, 0) 
magnetic excitations expected from our theoretical proposal. 
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SUPPLEMENTARY MATERIAL - 
ANTIFERROQUADRUPOLAR AND ISING-NEMATIC 
ORDERS OF A FRUSTRATED BILINEAR-BIQUADRATIC 
HEISENBERG MODEL AND IMPLICATIONS FOR THE 
MAGNETISM OF FESE 

The ground-state spin configurations in the classical spin model 

Exactly at T = 0, the classical 0(3) spins are always or¬ 
dered. Therefore, the AFQ order is accompanied by magnetic 
dipolar orders. Because the {n, 0) AFQ order doubles the unit 
cell, the structure factor Sd( q) of the compatible magnetic 
dipolar order must show a two -q structure as the consequence 
of Brillion zone folding, i.e., 5/j(q) = 5 n (q + q..\ ) where 
q .4 = ( 74 , 0 ). The ordering wavevector q depends on model 
parameters. 



(b) 


(d) 



FIG. SI. Four among the infinitely degenerate ground-state spin pat¬ 
terns in the case of J 3 = 0 and Ji = J 2 . In each case, the dashed box 
shows the magnetic unit cell. The corresponding ordering wavevec- 
tors are as follows: q = (0, n) and (n, ^ r) in (a); q = (± 7 t/ 2 , 7 r) in 
(b) and (c); q = (27r/3,7r) and (— 7 r/ 3 , n) in (d). 

We find that in the (tt, 0)/(0, 7 r) AFQ ground state the spins 
are ordered at a wavevector ( q, 7 r)/(7r, q) with infinite degen¬ 
eracies for J 3 = 0. Assuming a ( 7 r, 0) AFQ order, the spin 
variable at site i is S* = £( i x )[e x ^(i x ) + e y (l - l(ix))]l{iy), 
where i x and i y are coordinates of site i, r ){i x (y)) = (1 + 
(_l)Mjh)/ 2 , and £,(i x ) = ±1 is a random variable defined 
on each column of the lattice. The randomness in the real- 
space spin configuration leads to infinite number of degener¬ 
ate ground-state spin patterns. Transforming to the momen¬ 
tum space, they correspond to ordering wavevector at (q. 7r) 
(and (q + 7 r, n)) with q an arbitrary number. Some of the de¬ 


generate spin patterns are shown in Fig. SI 


As for the case of J\ = 0 and J 2 = J 3 , in the AFQ phase, 
we still find 16-fold degenerate ground-state spin patterns at 
ordering wavevectors (±7r/2, ±7 t/2 ). Some of the spin pat¬ 
terns are shown in Fig. [S2] In both cases, the large number 
of degenerate classical spin ground states helps to stabilize an 
AFQ order without a magnetic dipolar one when the quantum 
fluctuations are taken into account. 




t^J 


FIG. S2. Two out of the 16 degenerate ground-state spin patterns 
in the case of Ji = 0 and J 2 = J 3 . In each case, the spin pattern 
can be obtained by repeatedly aligning the spins in the dashed box 
in a staggered way. In both cases, the ordering wavevectors are q = 
(±7t/2,±7t/2). 


Another interesting observation is that at T = 0, the dom¬ 
inant ( 74 , 0) AFQ order of Q x ~ y coexists with a sublead- 
ing FQ order of Q r . This can be checked using the 
ground-state spin configurations in Figs. [ST] and |S2| which 
gives Q x ~ 3z = l/-\/3 at each site and Q x y2 = ±1. 


Spin excitations and Goldstone modes in the quantum 5=1 
model 


For the case of quantum spin 5 = 1, the model defined in 
Eq. (1) of the main text can be studied by an SU(3) Schwinger 
boson approach.(2 At each site, let | — 1), |0), and |1) be the 
three eigenstates of the spin operator S z . We can define a 
time-reversal invariant basis of the SU(3) representation: 


|Z> = ^(|1>-|-1»: 

|y) = -T(|l) + |-l) )! 

I z) = -i|0). 


(SI) 


Within this representation, we can then define three 
Schwinger bosons associated with the above three states, 
*4,10) = I a), where a = x,y, z, and 10 ) is the null state of the 
Schwinger bosons. The three bosons satisfy a local constraint 
at each site: 


E 6 


bjrv — 1 - 


(S2) 


The spin dipolar and quadrupolar operators can be written in 
terms of the Schwinger boson bilinears as 


(S3) 


Q? = ~(blb ip + b\ p b ia ), 

Q x .~ v = ~(b\ x b ix - b\ y b iy ), 


(S4) 


Q r :~ iz = - b\ x b ix b] y b iy ), 
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where a , /3, and 7 run over a:, y, and 2 , and e a ^ 7 is the Levi- 
Civita symbol. The Hamiltonian is then rewritten as 

H = 2 E/ {jnbl a bj a bjpbip + (K n — Jn)b\ a tfj a bjpbip^ , 
i,5 n ,Ot,(3 

(S5) 


where j = i + and (with n = 1, 2, 3) connects site i 
and its n’s nearest neighbor sites. 

We assume the following ground state at the mean level: 
\ip) = Ui\^i) = n, J2 a /»a&Ll 0 )> where the coefficients sat¬ 
isfy Yl a \ fia \ 2 = 1- The (77 0) AFQ order can be obtained 
by requiring condensation of b x and b y bosons at sites in odd 
and even columns, respectively. Correspondingly, the mean- 
field ground-state wave function at site i is \ipi) = \x)i if the 
x coordinate of site i ( i x ) is odd, and | t/ 7 ) = |j/) i if i x is 
even. One could check that this wave function is indeed as¬ 
sociated to an AFQ order at wavevector q,i = (77 0) since 

(V’il'EEV’i) = e qA ' ri = ±1. 

We study the spin excitations in the AFQ phase by using 
a flavor-wave theory. We first perform a local rotation in the 
spin space. 



cos 0i sin 0i 0 \ 
— sin &i cos 6i 0 
0 0 1 / 



(S 6 ) 


such that in the rotated basis, only one flavor of bosons, d ix , 
condenses. In the ( 7 r, 0) AFQ phase, this corresponds to tak¬ 
ing di = 0 if i x is odd and 9i = 7 t/2 if i x is even. Using the 
constraint in this rotated basis, d\ a di a = 1 , we obtain 


dix > y 1 Aydiy d\ z di Z 

~ 1 - \ (d\y d iy + d \z d iz)- (S7) 

Using Eq. ( |S7[ >, we can expand the Hamiltonian in Eq. ( |S5| > 
in terms of the magnon operators di y , d zz , and their Hermitian 
conjugates. We then truncate the expanded Hamiltonian to 
keep up to the quadratic terms of di y , di Z , and their Hermitian 
conjugates. Given that the ground state is the AFQ state, the 
linear terms in di y , di Z automatically cancel out, and we arrive 
at, up to a constant energy, a quadratic Hamiltonian. After 
Fourier transforming it to the momentum space, this quadratic 
Hamiltonian reads, 

H ~ ^ ' A ka (d ka dka T d—ka d — ka ) 

k,ac=y,z 

+Bka( d \. a d ' l _ ka + dkad-ka). (S 8 ) 


Here k runs over the unfolded Brilluion zone (BZ), and 

Ak y = J\ (cos k x + cos k y ) — Ki cos k x 

— 2(I\ 2 — J 2 ) cos k x cos k y + J 3 (cos 2 k x + cos 2 k y ) 

+ 2K 2 ~2K 3 , (S9) 

B ky = Ki cos k y — Ji(cos k x + cos k y ) — 2 J 2 cos k x cos k y 
+ (K 3 — J 3 )(cos2k x + cos2k y ), (S10) 

Akz = J\ cos k y + J 3 (cos 2 k x + cos 2 k y ) — I\\ — 2K 3 , 

(Sll) 


Bkz = (K 1 — Ji) cosky + ( K 3 — J 3 )(cos2k x +cos2 k y ). 

(S12) 


The Hamiltonian in Eq. ( |S 8 | > can be diagonalized via a Bo- 
goliubov transformation 

d ky{z) ~ u fcy(z)7fci/(z) ~ v ky(z)'Yky(z)’ (S13) 

and 

H « E £ fc«(7L7fea+7 — ka.'Y — hot ), fS 14) 

k,ct=y,z 

where the magnon dispersion 


and 


€kcx — 



B 2 


V'ka 


1 ( A ka 



i 

'Vkot. — 2 


A 


key 


VA. 


2 — B 2 

ka ^ key 



(S15) 


IS 16) 


The parameter regime where the (77 0) AFQ phase is stable 
is obtained by requiring A ka 7 0 and ,4/. — B ka ^ 0 at 
every k in the entire BZ and for both a = y,z, and can be 
determined numerically. 

The dynamical spin dipolar and quadrupolar structure fac¬ 
tors SQ 1 (q 1 uj) and S'jj Q (q, oj) can be calculated within the 
diagonalized representation. In gerneral. 


S^( q,w) = Y, \(9\QZH\ 2 6(u - (E m - E g )), (S17) 

m 

S% a ( q, w) = E \(9\SZ\m)\ 2 5(u; - (E m - E g )). (S18) 

m 

Here, | g) and \m) refer to the ground state (with eigenenergy 
Eg) and the to’s excited state (with eigenenergy E m ) in the 
flavor-wave theory. 


Qq 

QCt 

Oq 


Vn 


E- 


,iqr 






i 

^i > 


(S19) 


(S20) 
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where N is the number of spins of the system, and and Sf 
are expressed in terms of Schwinger bosons using Eqs. 
and (|S4|>. For example, for /r = x 2 — y 2 , 



2 

^2 d ky dk+q+q A ,y 

^lz^fc+q+qA,z- (S21) 


This leads to Sq 1 ( q, w) = N8q,q A 5(uj) up to the one-magnon 
contribution, which confirms the AFQ order at q A = (n, 0). 

The one-magnon contribution to the spin dipolar correlation 
function is also non-zero. We find that the transverse dynami¬ 
cal structure factor 

Sf?(q, W ) = S™(q,u>) = - £qz ) 

(S22) 

From Eqs. ( |S1 l| l,( [ST2| ), and ( |S15| ), we find a Goldstone mode 
near q, 4 = ( 7 r, 0). For q = q A + (dq x ,dq v ), 

£q Z « \Jvldql + Vydqy, (S23) 

with anisotropic velocities 

v x = 2y/K 3 r], 

v y = y/{K 1 + 4 /v 3 )t 7 , (S24) 


where r) = K] + 2 K 3 — ,/| - 2 J 3 . The Goldstone mode near 
q i is also seen in the dynamical spin dipolar structure factor: 

Sff(q. <Ia,u) ~ - Cq z ). (S25) 

Note that the static dipolar stmcture factor Sff (q A ) = 0, be¬ 
cause of the absence of long-range magnetic dipolar order in 
the AFQ phase. But the Goldstone modes associated with the 
broken spin rotational symmetry can be observed from the dy¬ 
namical spin dipolar structure factor. A representative plot of 
S'ff‘ (q, oj) showing the Goldstone modes is displayed in Fig. 4 
of the main text. 


Quantum fluctuations and coupling to itinerant fermions. 

The field theory that describes the two coupled order pa¬ 
rameters Q 4 and Q/ ; will be similar to that of the ( 77 , 0) AFM 
order of the J \ — J 2 Heisenberg model. The effect of coupling 
to the itinerant fermions can be treated as in Ref. 0 within an 
effective Ginzburg-Fandau action: 

S(Q a Mb) = J d{q} J d{u>}[S a +S* + ...], 

5 2 (q,w) = r(w,q)[|Q A (q,w)| 2 + |Q s (q,w)| 2 ] 

+ v(q 2 x - g 2 )[Q A (q,w) • Q s (-q, -w)], 
s 4 ({q}, M) = «(|Qa| 4 + |Qs| 4 ) + m'|Q a | 2 IQs| 2 

+ w|Qa-Qs| 2 . (S 26 ) 


where r(w, q) = Vq + wA + w 2 + 7 |w| + cq 2 , r 0 < 0 and 

A > 0 are constants, w is the coherent quasiparticle spectral 
weight of itinerant electrons, and 7 is a Fandau-damping co¬ 
efficient. Note that u < 0, and {q}, {u;} mark the set of four 
momenta and four frequencies that enter S ,\; the momentum 
and frequency integrals are understood to each contain a delta 
function that fixes the sum of the momenta and the sum of the 
frequencies at zero. A similar form also exists for the AFM 
orders, M A and M 0 . The shift of r by w and the damp¬ 
ing may lead to the loss of magnetic order in the system 0 , 
and stabilize a pure AFQ order. In the absence of the AFQ or¬ 
der, the Ising-nematic order will be concurrently suppressed 
0 . However, when the dominant order is AFQ, one can read¬ 
ily reach the regime where the quantum fluctuations eliminate 
the weaker AFM order while retaining the stronger AFQ order 
and the associated Ising-nematic order. 


The evolution of the Ising-nematic order parameter and 
transition temperature as a function of J2/A2 


Here we study how the dominant Ising-nematic order pa¬ 
rameter and the associated transition temperature, T a , change 
with varying the J 2 IK 2 ratio in the 2D classical bilinear- 
biquadratic Heisenberg model. The calculations are done via 
classical Monte Carlo simulations using Metropolis sampling 
with up to 80 x 80 lattices and 10 5 Monte Carlo steps. 

We tune the J 2 IK 2 ratio by going along the blue dashed 
trajectory in the phase diagram shown in Fig. |S3fa). As dis¬ 
cussed in the main text, with increasing K 2 . the ground state 
of the system changes from the ( 7 r, 0)/(0,7r) AFM to the 
( 7 r, 0) /(0,7r) AFQ state. We have tracked the evolution of the 
Ising-nematic order parameters 04 and 04 , and find that the 
change of the ground state is reflected in the variation of the 
Ising-nematic order parameters: the dominant Ising-nematic 
order parameter changes from 04 in the AFM phase to <r 2 in 
the AFQ phase, as clearly shown in Fig. |S3f b) and (c). We 
have also determined the transition temperature T a from our 
numerical results. T a first decreases then increases with in¬ 
creasing K 2 along the trajectory, showing a remarkable min¬ 
imum near the phase boundary between the AFM and AFQ 
phases (Fig. |S3[d)). This minimum indicates enhanced fluc¬ 
tuations around the Ising-nematic order due to the competition 
between the AFM and AFQ orders. Our results on the evolu¬ 
tion of the dominant Ising-nematic order parameter and T a 
with J 2 / K ‘2 are fully consistent with the general picture pro¬ 
posed in Fig. 5 of the main test. When an interlayer coupling 
is turned on, we expect similar results for the Ising-nematic 
order, because the change of the dominant Ising-nematic or¬ 
der parameter and the evolution of T a reflect the competition 
between the underlying AFM and AFQ ground states. In this 
case, in the dominating AFM regime, it is well known that a 
nonzero transition temperature develops for the AFM order. 
Similar reasoning applies to the dominating AFQ regime. 
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K 2 / l K ll 


FIG. S3, (a): Mean-field ground-state phase diagram of the clas¬ 

sical bilinear-biquadratic Heisenberg model at J 3 = K 3 = 0 and 
Ji = J 2 . The blue dashed line shows the trajectory by tuning 
J 2 /K 2 , and the red and black dots show two representative points 
along this trajectory, (b) and (c): Temperature dependence of the 
Ising-nematic order parameters g\ and a 2 at the two representative 
points of the phase diagram in panel (a). There is a change of the 
dominant Ising-nematic order parameter from <j\ to 02 when the 
ground state varies from the AFM to the AFQ phase. The blue dashed 
line marks the transition temperature, T„, of the Ising-nematic transi¬ 
tion. Data shown are based on Monte Carlo simulations on a 40 x 40 
lattice, (d): Evolution of T a with varying J 2 / K 2 along the trajec¬ 
tory shown in panel (a). A minimum of T a is located near the phase 
boundary (black solid line) between the AFM and the AFQ phase, 
which is consistent with the proposal made in the main text, as illus¬ 
trated in Fig. 5 of the main text. 
















